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Quantum transport of strongly correlated fermions is of central interest in condensed matter 
physics. Here, we present first-principle nonequilibrium Green functions results using T-matrix 
selfenergies for finite Hubbard clusters of dimension 1,2,3. We compute the expansion dynamics 
following a potential quench and predict its dependence on the interaction strength and particle 
number. We discover a universal scaling, allowing an extrapolation to infinite-size systems, which 
shows excellent agreement with recent cold atom diffusion experiments [Schneider et al., Nat. Phys. 
8, 213 (2012)]. 
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Particle, momentum and energy transport of strongly 
correlated quantum systems are of growing current in¬ 
terest in condensed matter [1-3], ultracold quantum 
gases [4-6] and dense plasmas [7]. Recently, direct mea¬ 
surements of quantum transport based on the expan¬ 
sion of ultracold atoms following a confinement quench in 
Hubbard-type one- and two-dimensional (ID, 2D) opti¬ 
cal lattices with single-site resolution were reported [4, 6]. 
Also, the dynamics following a quench in lattice depth 
have been measured [8-11]. On the other hand, theoreti¬ 
cal studies of these transport processes pose fundamental 
challenges. The authors of Ref. [4] presented numerical 
results for the expansion of ultracold fermions in a 2D 
lattice from a semi-classical Boltzmann equation model 
with a collision integral in relaxation time approxima¬ 
tion (RTA) and reported good overall agreement with 
the experiment. At the same time they pointed out that, 
while the expansion can be modeled in ID using a density 
matrix renormalization group (DMRG) approach [5, 6], 
“[...] so far no methods are available to calculate the 
dynamics quantum-mechanically in higher dimensions.” 

It is the purpose of this Letter to fill this gap [12]. 
We present first-principle nonequilibrium Green func¬ 
tions simulations using the full T-matrix approxima¬ 
tion applied to a fermionic Hubbard model of dimension 
D = 1... 3. We analyze the expansion of an initially 
conhned system and resolve the short-time dynamics of 
the particles demonstrating how correlations build up in 
finite strongly correlated inhomogeneous fermionic sys¬ 
tems. We also analyze the long-time limit of the expan¬ 
sion velocity, u“p, for ID, 2D and 3D systems, for a 
broad range of coupling strengths U and particle num¬ 
bers N and observe a universal scaling u“p ~ 1/Vn, 
independent of U and D. This enables us to extrapolate 
our results to the macroscopic limit and directly com¬ 
pare with the measurements of Ref. [4]. The agreement 
is excellent and achieved without any free parameters. 


Nonequilibrium Green functions (NEGF) in T- 
matrix approximation (TMA). The NEGF are de¬ 
fined on the complex Keldysh contour C with time¬ 
ordering operator Tc as 

, ( 1 ) 

for lattice site indices s = {si,..., sd), s' and spin 
projection cr G {t)'!-}- The equations of motion for 
the NEGF are the Keldysh-Kadanoff-Baym equations 
(KBE) [13, 26], 

( 2 ) 

= Sc{z- z')'^s,s' + dzE^-(z, z)Gl^,{z, z '), 

and its adjoint (summation over s is implied). The corre¬ 
lation part of the selfenergy (in excess to Hartree-Fock) 
is computed on the T-matrix level and reads, for the 
Hubbard model [16], 

z') = z') Gt[p(z', z), (3) 

T^Az,z') = -ihU^Gl,{z,z') (4) 

+ ihU J^dzGl-iz,z)T,^,iz,z'), 

Gl,iz,z')=GlAz,z')Gt,A^')- 

Here, T can be understood as an effective interaction 
obeying the Lippmann-Schwinger equation (4), e.g. [IS¬ 
IS]. Taking only the leading term in (4), which describes 
the interaction with a single electron pair, the TMA re¬ 
duces to the second-order Born approximation (SBA). 
We underline the conserving character of this approxi¬ 
mation [13] and, in fact, conservation of particle number 
and total energy is observed to high accuracy in all our 
simulations [27]. 
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Figure 1. (Color online) Time evolution of initially circularly 
confined fermions in a two-dimensional 19 x 19 Hubbard lat¬ 
tice for the three particle numbers N = 2, 26, 72 for U = 1 
(upper three rows) and [7 = 4 (bottom row). Results of TMA 
calculations. Columns correspond to four representative time 
steps. Color code corresponds to 


The use of this complex approximation under full 
nonequilibrium conditions has only recently become pos¬ 
sible for the Hubbard model, e.g. [16-18]. Compar¬ 
isons with exact diagonalization calculations (Cl) con¬ 
firmed the high accuracy of this approximation [17, 19]. 
We performed additional simulations using T-matrix 
and second-order Born selfenergies with the general¬ 
ized Kadanoff-Baym ansatz (GKBA) with Hartree-Fock 
propagators to reduce selfconsistency effects that are 
known to be critical in finite systems [16]. 

Numerical results. The two-time KBE (2) were 
solved for a H-dimensional Hubbard model of sites 
and a step-like circular confinement V^{t) of radius R 
(switched-off for t > 0, as in the experiment [4]), 


H{t) 


{s,s') s 

+ J2J2 ( 5 ) 

s 


where (s, s') denotes nearest neighbor sites. It was 
studied for H = 1... 3 and a broad range of coupling 
parameters, 0 < U <8. The number of fermions, 
N = + N^ = 2N^, was varied in the range 2... 114. 

The calculations started in the thermodynamic ground 
state which agrees well with the experimental conditions 
[4]. After the confinement was switched off (potential 
quench), the diffusion of the fermion cloud was recorded. 
Figure 1 shows snapshots of a typical expansion for two 
couplings, [/ = 1, 4. For U = 1, the density rapidly 


evolves towards square symmetry of the lattice, whereas 
for {7 = 4 the core region remains circular, as observed 
in the experiment [4]. It is apparent that, there is first 
a universal initial phase where diffusion is only possible 
for particles at the cluster edge due to Pauli blocking 
[5, 17, 20]. This is followed by a more complex evolution 
that strongly depends, both, on U and N (compare the 
rows in Fig. 1). 

To quantify this evolution, we introduce the cloud di¬ 
ameter d, corrected for the initial cloud diameter, i?^(0), 

dit) = v/i?2(t)-i?2(0), (6) 

s s 

that involves the time-dependent site occupation num¬ 
bers n^(t) and the static cloud center-of-mass, Sq. The 
upper part of Fig. 2 shows the dynamics of the instan¬ 
taneous expansion velocity, Vexp{t) = for various 

U, for the test case N = D = 2 where exact diago¬ 
nalization (Cl) data are available. Our NEGF results 
within TMA and GKBA-l-TMA show good agreement 
with Cl [28] which gives us confidence in their reliability 
for much larger systems and higher dimensions that are 
out of reach for exact methods. 

The expansion velocity starts from the ideal ballistic 
value, 'yexp(O) = V2D = 2, and converges to a constant 
smaller value, that monotonically decreases with U. 
The dynamics of Vexp{t) between these limits are result¬ 
ing from the non-trivial interplay between single-particle 
and correlation effects which are straightforwardly acces¬ 
sible within our NEGF approach. Of particular interest 
are the single-particle and correlation energy {E^p, Ecorr) 
[17] as well as the entanglement entropy S = Ssp + Scorr 
[ 21 , 22 ], 



-(l-n^+ log 2 (l - n, + > (7) 


where is the double occupation 

of site s. The single-particle part, Sgp, follows from the 
replacement n^J), —)■ in Eq. (7), and the correlation 

part is the remainder. The dynamics of the energy and 
entropy contributions allow us to identify three charac¬ 
teristic phases of the evolution: during the first, ^sp and 
Esp are built up, resulting in a decrease of t^exp, see top 
part of Fig. 2. Here, the increase of Ssp measures the 
transition from a state of independent particles {S = 0) 
to an interacting many-body state. The inflection point 
Tgp (circles) of Ssp (and Esp) is representative for the 
time scale of this phase. Subsequently, a second phase 
ensues that is characterized by the saturation of Esp lead¬ 
ing to a convergence of Vexp- The simultaneous build-up 
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Figure 2. (Color) Top left: Time evolution of the expansion 
velocity for the N — 2 setup of Fig. 1 and five values of U. 
The horizontal extrapolation of Uexp is used below to charac¬ 
terize the stationary expansion of the system. Different line 
styles denote the TMA, GKBA-I-TMA and the exact result 
(Cl). Symbols: inflection points of the characteristic evolution 
quantities (see top right figure for explanation). Top right: 
Exact evolution of the single-particle and correlation parts of 
the entropy and energy, for N = 2 and ?7 = 4. Symbols mark 
the respective inflection points (positions agree with TMA 
results (not shown here)). Shaded areas correspond to the 
three phases of the evolution (see text). Bottom: Inflection 
times Tsp (single-particle, solid lines with diamonds) and Tcorr 
(correlation, dashed lines with circles) from TMA. right: De¬ 
pendence on U for N = 2,58,114. left: Dependence on N 
for U = 1.0, 2.5,4.0. 

of correlations partly prolongs the saturation and deter¬ 
mines the final value of Uexp- The representative time 
for these processes, Tcorr> is the inflection point of S'corr 
(and Ecori, diamonds). In the third phase, the expan¬ 
sion velocity and correlations are saturated, whereas the 
single-particle observables continue to increase with the 
ongoing expansion. 

Both characteristic time scales show an interesting de¬ 
pendence on U, cf. upper left part of Fig. 2, and also on 
N. This is further explored in the bottom part of Fig. 2, 
where we show Tgp and Tcorr for varying U (left) and N 
(right). It is evident that Tcorr is one order of magnitude 
larger than Tsp—in striking contrast to homogeneous sys¬ 
tems [23], and both increase with N. The reason is that 
entanglement entropy (and energy) are first produced at 
the cluster edges and the build-up continues towards the 
center once the outer doubly occupied sites are depop¬ 
ulated. For small N, the active regions quickly overlap 
impeding the production. On the other hand, when U is 
increased, the effective scattering rates increase what ac¬ 
celerates the inward propagation of energy and entropy. 

To quantify the stationary properties, we extrapolate 
the expansion velocity to t —>■ oo. The result v'^^{U]N) 



Figure 3. (Color online) Left: Asymptotic expansion velocity 
vs. particle number for D = 1... 3 and U = 1... 3. Symbols 
correspond to actual TMA data points, errors are smaller 
than symbol size. Dashed lines: linear extrapolation N —^ 
oo according to Eq. 8. Right: Corresponding slope y vs. 
bandwidth-normalized interaction. 


is well suited to characterize the influence of correlations 
and of finite-particle number effects. Let us now consider 
how the diffusion properties depend on the cloud size. To 
this end, we analyze the asymptotic expansion velocity, 
v^p{U; N), at fixed U, for varying N. The details of 
the extrapolation and the calculation of the error are 
explained in the supplementary material [30]. The left- 
hand part of Fig. 3 shows the obtained results for U = 
1... 3 and different dimensions D = 1... 3. Beside the 
above-mentioned decrease of with U, it also shows 
an increase with D, that is due to the enlarged number 
of degrees of freedom. A striking observation is that, for 
all U and D, the asymptotic expansion velocity exhibits 
a linear scaling with 

Cp(C/; N; D) - D) = x(C/; D)N-^/^ , (8) 

for sufficiently large N. The right-hand part of Fig. 3 
shows the dependency of the slope x on the bandwidth- 
normalized interaction strength C//(6/2) with the effec¬ 
tive bandwidth b = AD. For all dimensions, x starts 
at zero for vanishing U, increases to a maximum below 
U = (6/2) and decreases again for further increased U. 
The A^-independent value x(0, D) = 0 is a consequence 
of ballistic expansion of independent particles. On the 
other hand, for {7 —?> oo, the doubly occupied sites are 
effectively frozen and the particles do not expand regard¬ 
less of N. In-between these two limits the slope shows a 
qualitatively similar behavior for all dimensions, exhibit¬ 
ing a steep rise, for small U and a slow decrease, for large 
U. 

To further explore the universality of the scaling with 
we transform the Hamiltonian to momentum 
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Figure 4. (Color online) Momentum distribution p{k) at 
t — 9.5 for a ID system of Ns = 65 for [7 = 3 and differ¬ 
ent particle numbers N = 2 ... 42 obtained with TMA. The 
dashed line denotes the initial distribution po{k) which is uni¬ 
form and independent of N. Inset shows the amplitudes a of 
the distribution vs. particle number for U = 2 ... 4. Symbols 
correspond to actual data points, lines denote linear hts. 

space [24], 

+ ^ ^p+q,t^k-q,l^p,t^k,l^ (^) 

^ k.p.q 

where e{k) = —2 cos{ki) and V^{t) is the transform 
of Figure 4 shows the momentum occupation 

probability p{k) = n{k)/N of a ID system for f/ = 3 
and 2 < < 42 at the end of the simulation. For all TV, 

p{k) oscillates around a constant mean with an ampli¬ 
tude a, that monotonically decreases with TV. For large 
TV, we observe 

p(k) =p(fc) = ^-bacos(A:), (10) 

where the value of a({7, N) is shown in the inset of 
Fig. 4 for different U and we again encounter a scal¬ 
ing a ^ The recovery of the same asymptotic 

scaling for p{k) as for the expansion velocity is another 
indication of universal behavior. This coincidence is not 
surprising since p{k) directly enters the kinetic energy 
and also determines the interaction energy. 

The observed robust scaling makes us conhdent to 
use the extrapolation of v^p{U;N) with respect to N 
for quantitative predictions of the diffusion properties of 
strongly correlated macroscopic systems. In the inset of 
Fig. 5 we show the result 14xp(b^) = limAr_).oo 
as a function of U and confirm the monotonic reduction 
of the mobility with the interaction strength. The HF 
approximation exhibits strong deviations which under¬ 
lines the key role of correlations. Thus, having obtained 



Figure 5. (Color) Extrapolated core expansion velocity Cexp 
from NEGF (circles with error bars), compared to the ex¬ 
periments (plus-signs) for different recoil energies Er [29] and 
the relaxation time approximation results (gray dashes) of 
Schneider et al. [4]. Inset: NEGF results for the asymptotic 
expansion velocity Fexp—TMA vs. HF and SBA. 


macroscopic results allows us to compare with the ex¬ 
perimental data that are typically obtained for large sys¬ 
tems with N ~ 10® fermionic atoms. To this end, we 
compute the core expansion velocity, Cexp, introduced in 
Ref. [4]—the velocity of the half width at half maximum 
(HWHM). In the experiments, Cexp showed an interest¬ 
ing zero crossing, around C = 3, after which it became 
slightly negative, indicating shrinkage of the central part. 
Even though the “core” may be an ambiguous definition, 
it is useful for quantitative comparisons. In fact, our 
T-matrix NEGF results show a surprisingly good agree¬ 
ment with the experimental data in the entire [/-range, 
including the zero crossing at [/ « 3, cf. Fig. 5. 

To summarize, we have presented first-principle NEGF 
results for finite ensembles of strongly correlated fermions 
on ID, 2D as well as 3D Hubbard lattices. With the 
developed simulations within the T-matrix approxima¬ 
tion, the expansion dynamics following a confinement 
quench could be accurately simulated. These dynamics 
were quantified using the expansion velocity and the rel¬ 
evant energy and entropy contributions, by which their 
coupling and particle number dependence were revealed. 
Our results show that full two-time quantum simulations 
can be successfully applied to mesoscopic fermionic sys¬ 
tems and, via extrapolation, even can be extended to 
macroscopic systems. The agreement with experimental 
results is excellent taking into account that our theory 
has no free parameters. Our approach is directly applica¬ 
ble to other transport quantities including electrical and 
heat conductivity and magnetic properties. Furthermore, 
the fully inhomogeneous character of the simulations al¬ 
lows one to study the influence of the system geometry 
and dimension. Finally, it will be interesting to experi¬ 
mentally verify the observed particle number dependence 
of the expansion velocity and other transport properties. 
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EXTRAPOLATION OF THE ASYMPTOTIC 
EXPANSION VELOCITY 


The asymptotic expansion velocity w“p is obtained 
from Vexp{t) by averaging after reaching convergence at 
the time tavg) which fulfills 



dUexp 

dt 


it) 


< e, 


( 1 ) 


for all t > tavg and a given small parameter (e <C 1). To 
quantify the error of u“p, we use the standard deviation 
a (f^p) of the averaging process. The macroscopic limit 
Hexp is extrapolated via 


^;Sp(iV) = Vexp + xfV-'/^ (2) 


where T4xp and the slope y are fit parameters. Only 
particle numbers N larger than a cutoff value N* are 
taken into account. The errors a (r’^p) are also included 
in the fit process, resulting in the final accuracy a (T4xp)- 


COMPUTATION OF THE CORE EXPANSION 
VELOCITY 

The core expansion velocity c“p and its macroscopic 
limit Cexp are obtained similar to [1]. For D = 2, the 
density distribution is averaged azimuthally for each time 
step. The half width at half maximum i?HWHM(^) of the 
resulting profile is used to measure the core width. Ad¬ 
justed for the initial core width, o^p deter¬ 

mined by fitting the resulting i?HWHM(0 to (cf. Eq. (6) 
of the main text) 

Rnwnuit) = \J (A^whm)^ + (cSpt)^ (3) 

for all t > tavg with i?HWHM ‘^Sp fitting pa¬ 

rameters. Since the core of the density distribution starts 
to shrink for sufficiently large interaction strength U, we 
apply 

.RHWHM(t) = \/(^HWHm) ~ (cSpt) (4) 

instead, following [1], and consider c“p the speed of con¬ 
traction of the core region. 


The extrapolation of c“p is done as in Eq. (2), re¬ 
sulting in the macroscopic core expansion velocity Cexp, 
confirming the scaling with 1/^/IN). 


DEPENDENCE OF THE EXPANSION 
VELOCITY ON THE SHAPE OF THE INITIAL 
CONFINEMENT 

One may wonder whether the results of the main paper 
depend on the chosen steep confinement, cf. Fig. 1 of the 
main text. Here, we demonstrate that the shape of the 
initial confinement has only a minor effect on the asymp¬ 
totic expansion velocity. To this end we use a harmonic 
confinement 


V{R) = , (5) 

with curvature ■jk ■ To achieve a similar shape for different 
N, we choose 


jk{N) = k/N, (6) 

for three strengths fc = 3, 5,10. Together with the step¬ 
like potential (fc —)■ oo), the results are shown in Fig. 1. 
Even though the initial density profile is affected by the 
curvature (see inset), the expansion velocity is not. In 
particular, the macroscopic limit, Uexp changes by less 
than 10%. 


COMPARISON OF DIFFERENT MANY-BODY 
APPROXIMATIONS 

To assess the influence of the chosen selfenergy approx¬ 
imation on the expansion, we computed the expansion 
velocity for a representative N = 58, D = 2, U = 2.5 sys¬ 
tem, shown in Fig. 2. While, in Hartree-Fock approxima¬ 
tion, Vexp only slightly decreases compared to the ideal 
velocity, the correlated methods accordantly capture a 
steeper decline and a later inset of convergence. Ana¬ 
lyzing the initial behavior, one notices that, at first, the 
TMA and SBA start to deviate from each other. Fur¬ 
ther on, each branch divides with respect to the GKBA. 
In SBA, GKBA-I-SBA and TMA, Vexp converges to a 



2 



N-1/2 


value which differs by less than 10%. In comparison, 
the GKBA+TMA show a larger deviation. To classify 
the reliability of the approximations, we state that the 
TMA is superior to SBA by construction [2]. Regarding 
TMA and GKBA+TMA, which are both strong coupling 
approximations, we have shown that they differ only in 
the degree of selfconsistency [3]. The N = 2 test against 
exact results (cf. main text, Fig. 2) and the agreement 
with the experiment (cf. main text, Fig. 5) slightly favor 
TMA over GKBA+TMA, though. However, no ultimate 
conclusion on the preferable method can be drawn for 
the considered systems and quantities. 


Figure 1: Dependence of the asymptotic expansion velocity 
on N for confinement potentials of different curvature 7 fe(A). 
Insets show the shape of the initial density profile. 
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Figure 2: Time dependence of the expansion velocity for a 
typical example of A = 58, D = 2 and U = 2.5. Compari¬ 
son of selfenergy approximations: Hartree-Fock(HF), second- 
order Born (SBA), T-matrix (TMA), with and without addi¬ 
tional generalized Kadanoff-Baym ansatz (GKBA). 






















































